Universality in three-dimensional Ising spin glasses: Nonequilibrium dynamics 

from Monte Carlo simulations 
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The non-equilibrium dynamics of the three-dimensional Edwards-Anderson spin-glass model with different 
bond distributions is investigated by means of Monte Carlo simulation. A numerical method is used to determine 
the critical temperature and the scaling exponents of the correlation and the integrated response functions. The 
results obtained agree with those calculated in equilibrium simulations and suggest that the universality class 
does not depend on the exact form of the bond distribution. 
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The critical behavior of spin glasses is a fundamental sub- 
ject of interest in the statistical mechanics of disordered and 
frustrated systems. The real nature of the phase transition 
is still under discussion and the controversy about different 
related topics has not been resolved satisfactorily. In this 
context, it has been suggested that the basic universality rule 
which state that the critical exponents depend on the dimen- 
sion of space, the number of order parameter components and 
the range of interactions, does not hold in spin glasses.'"^ In 
particular, by means of a technique that combines equilibrium 
and dynamic measurements, it was shown that the universality 
class of the three-dimensional (3D) Edwards-Anderson spin- 
glass model^, depends on the exact form of the interaction 
distribution function. On the other hand, for several choices 
of bond distributions, the finite-size scaling analysis of the 
Binder cumulant, the correlation length and susceptibility cal- 
culated in equilibrium simulations, suggests that this model 
obeys universality.^ In addition, there is not satisfactory agree- 
ment between the critical temperatures Tc calculated by non- 
equilibrium and equilibrium techniques. 

In this work we propose a non-equilibrium method to deter- 
mine the r^ and the scaling exponents of the correlation and 
the response functions. We use this technique to study the 
3D Edwards-Anderson spin-glass model with three different 
bond distributions. The values the T^ that we obtain agree with 
those calculated in equilibrium simulations. Also, we found 
that the critical exponents are very close to each other, sug- 
gesting that the most probable scenario is that the universality 
class does not depend on the bond distribution form. 

The simplest non-equilibrium method to determine the crit- 
ical point is based on the temporal relaxation of the order 
parameter.^ First, at time f = the system is prepared in a fully 
ordered state and the dynamics is simulated with a standard 
Monte Carlo algorithm. Tc is estimated as the temperature at 
which, in the asymptotic regime, the order parameter follows 
a power law. This method is appropriate to study a wide vari- 
ety of systems. However, although the slow dynamics present 
in disordered and frustrated systems favors the application of 
non-equilibrium techniques, in general for spin glasses these 
do not allow an accurate determination of T^. The major prob- 
lem is that in simulations, different quantities seem to decay 
by a power law in a relatively wide interval of temperatures. 



To overcome this difficulty, additional scaling analysis of or- 
der parameter or susceptibility have been used to improve the 
resolution of these methods. ^'^'^ In addition, recently it has 
been proposed a different procedure based on the divergence 
of the relaxation time approaching the critical point. "^ 

Now, we will describe the non-equilibrium method pro- 
posed in this work. A typical protocol is used, which consists 
on a quench at f = from a disordered state (T ^ oo) to a 
low temperature T. From this initial condition we simulate 
the system with a Model A dynamics.'' Then, different two- 
time quantities are calculated which depend on both the wait- 
ing time t„, when the measurement begins, and a given time 
t > f,|,. In particular, for a system formed by A^ Ising spins, we 
determine the two-time autocorrelation function defined as 



Cit,tw) = - 



(=1 



(1) 



where (7; — ±1, (...) indicates an average over different ther- 
mal histories (different initial configurations and realizations 
of the thermal noise), and [...] represents an average over dif- 
ferent disordered samples. The scaling relation for C is 



C(f,fv.) 



'Mt/t^-) 



(2) 



where b is a non-equilibrium exponent.^-'— 

On the other hand, the corresponding two-time linear au- 
toresponse function is R{t,t„) ~ ^[^^j 5(a,(f))/5/;,(f,i-)], 
where /i,(fn.) is a time-dependent conjugated external field. 
We calculate in simulation the quantity 



p{t,t„) = r / duR{t,u), 
Jo 



(3) 



which is the integrated response when switching on the per- 
turbation only for times t < t„. The scaling relation for p is 



P (?,?«•) 



-7p(f /?«•), 



(4) 



where a is another non-equilibrium exponent. For critical sys- 
tems we have that b = a and, for t^. <C f — 1„, the scaling func- 
tions should follow a power-law decay, i. e., fc ^ [t jUv) '^'^ 
and /p ~ (f/f,v)^ '■'''', where A^. is the autocorrelation expo- 
nent and Zc is the dynamical critical exponent4ii^ 



The method proposed in this work consists in to calculate C 
and p for different values of T and t„ and, by means of a data 
collapse analysis, to determine the exponents b and a. Then, 
Tc is identified as the temperature for which the condition b = 
a is satisfied. As we shall see, this strategy will allow us to 
carry out a precise determination of the critical point. 

First, in order to validate the method, we study the two- 
dimensional (2D) ferromagnetic Ising model on the square lat- 
tice, for which Tc = 2/ln(l + \/2). The Hamiltonian is H = 
— Y.(i.j)J'^i'^j^ where the sum run over the nearest-neighbor 
pairs and J —1. A large lattice of linear size L = 300 (N ~L^) 
with fully periodic boundary conditions was simulated using a 
standard Glauber dynamics and, for each temperature studied, 
the averages were calculated over 5000 independent thermal 
histories. The correlation was calculated as usually, but the 
integrated response to an infinitesimal magnetic field was de- 
termined using the algorithm proposed in Ref. 13. This is very 
important to obtain a reliable value of the exponent a, making 
possible the realization of the method studied here. 

Figures[T](a) and[T](b) show, respectively, the data collapse 
of the correlation and the integrated response for different t„ 
and temperature T = 2.2692 sa T^, plotted as function of the 
variable x = f /f„, — 1 (we have used x instead of r/f„, but this 
choice has no consequence on the asymptotic behavior). The 
best data collapse was obtained by minimizing the sum of 
squared differences between all pairs of curves within a given 
range of x. For the range x > xq = 1, the exponents calculated 
in this way are b = 0.118 and a = 0.116. On other hand, if 
a common exponent c — b — a it is used to collapse simulta- 
neously both sets of curves, we obtain a value of c = 0.117. 
These numbers are very close to expected value for the Ising 
model, b = {d-2 + r])/zc ~ 0.115, where <i = 2 is the di- 
mension of the space, Tj = 1 /4 is the well-known static crit- 
ical exponent associated to the pair correlation function and 
Zc = 2. 1667. '2'^ As shown in Fig.[T](c), repeating this proce- 
dure for others temperatures the values of b and a cross at Tc, 
where an optimal collapse is obtained [see Fig.[T](d), where 
A^ is equal to the sum of the squared differences between all 
pairs of curves]. 

The previous example shows that the critical point can be 
located by looking for the temperature at which the condition 
b = a is fulfilled. This temperature can be regarded as Tc and 
the c value as a reasonable estimate of the non-equilibrium 
critical exponent b (or a). However, notice that for tempera- 
tures above or below Tc, only pseudo-exponents are obtained 
with this simple method. Close to (but not at) Tc, to calculate 
(when possible) the true non-equilibrium exponents of a given 
system, long-time simulations are necessary and maybe, ei- 
ther another scaling relations or a separation of the correlation 
and the response functions in their corresponding stationary 
and aging terms are required.'^) 

We now consider the 3D Edwards-Anderson spin-glass 
model whose Hamiltonian is, H = —Y.{i,j)Jii<^i<^i^ where 
the sum runs over the nearest neighbors of a cubic lat- 
tice with fully periodic boundary conditions. The bonds 
Jifs are independent random variables drawn from a given 
distribution P{Ji,j) with mean zero and variance one. As 
Ref. |4j, we concentrate on three distributions: the bi- 
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FIG. 1: (Color online) The 2D ferromagnetic Ising model. Panels (a) 
and (b) show, respectively, the data collapse of the correlation and 
integrated response at T = 2.2692 for different t„ as indicated, (c) 
The values of non-equilibrium exponents b and a, and (d) the quality 
of the collapses for different temperatures. 



modal PB{Ji.j) = [5{Jij-i) + 8{Ji.j+l)]/2, the Gaussian 
PciJij) — exp(y?^/2)/\/27r, and the Laplacian PiiJij) = 
exp(— \/2|7,-.y|)/\/2 distribution. In order to avoid confusions 
we will denominate, respectively, EAB, BAG and EAL mod- 
els to each one of these versions of the Edwards-Anderson 
model. Lattices of linear size L = 50(N = L^) were simulated 
using a standard single-spin Glauber dynamics and six values 
of f,, = 50, 100, 200, 400, 800, and 1600 were used. The in- 
tegrated response function was calculated as before for an in- 
finitesimal external field. The disorder average was performed 
over 3000 to 5000 different samples for each temperature. Be- 
cause it is expected that the scaling relations, Eqs. dU and (|4|i 
are valid for large values of the waiting time, we have studied 
different range of x and systematically we have discarded the 
curves with smallest t„. 

Figure 12] shows the results for the 3D EAB model. In panel 
(a) we can see the exponents b and a obtained by the present 
method for different ranges of x > xq, where all curves of C 
and p with t„ between 50 and 1600 were considered. Because 
the data agree reasonably well each other for xq > 3, from 
now on we will use xq = 3. Notice that the condition b — a 
is not fulfilled at any temperature. However, after discarding 
the data for f„, = 50 and next for f„, = 100, the Figs. |2](b) 
and[2](c) show that the curves begin to come closer together 
Finally, in panel (d) for the last three t„, we observe that 
the condition b = ais approximately fulfilled for T ~ 1.135, 
where we obtain b = 0.042, a = 0.046 and c = 0.044 (let us 
notice that it is not necessary that the curves cross to iden- 
tify the critical point). As this behavior is observed among 
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FIG. 2: (Color online) Non-equilibrium exponents as function of 
temperature for the 3D EAB model, obtained for different ranges 
of Ik as indicated. Curves in panel (a) were calculated for different 
values of .yq (arrows indicate how the curves change with increasing 
X()), while in (b), (c) and (d) jcq = 3 was chosen. 



T = 1.13 and T = 1.14, we conclude that T,. = 1.135(5) for 
the 3D EAB model. This value is very close to those obtained 
in the equilibrium simulations, e. g., Tc — 1.120(4) (Ref. 6) 
and Tc = 1.109(10),'*^ but is slightly lower than T,- = 1.17(4) 
(Ref. 9) and Tc — 1.19(1),''^ two critical temperatures calcu- 
lated from non-equilibrium simulations. Nevertheless, notice 
in Fig.|2]the tendency of the curves to merge atT — 1.17. Al- 
though a great number of samples were calculated, it was not 
possible to show that the condition b — a can be fulfilled at 
this temperature. In addition. Fig. [3] shows the data collapse 
of the correlation and the integrated response functions, where 
we have used b ^ a = 0.044. 

We have used the same protocol to study the others spin- 
glass models. Figure |4] (a) shows that for the 3D BAG 
model the condition b = a is approximately satisfied at Tc = 
0.95(1), where we determine that b = 0.046, a = 0.0455 and 
c — 0.0455. This condition is accomplished for ?„■ > 200. 
Again, this critical temperature agree very well with the value 
Tc = 0.95 1 (9) obtained in equilibrium simulations^ but, in this 
case, is slightly higher than Tc = 0.92(1), the value reported 
inRef.i4j 

On the other hand, Fig.|4](b) shows that for the EAL model 
is needed to discard the first four f„ to obtain a reliable value 
of Tf = 0.815(5). At this temperature we determine that 
b = 0.052, a ^ 0.042, and c = 0.047. The difference between 
these exponents is larger than the measure in previous mod- 
els. Probably, this is due to the fact that the Tc for this model 
is very low and bigger values of t and t„ need to be reached. 
Nevertheless, Fig.|4](b) shows strong evidence that this value 
corresponds to the true critical temperature which, as before. 
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FIG. 3: (Color online) Data collapsing of (a) the correlation and 
(b) the integrated response functions, for the 3D EAB model at 
T = 1 . 135. The value b = a = 0.044 was used. 



is found to be slightly higher than Tc ~ 0.72(2), the value re- 
ported previously.^ 

For comparison, we have also calculated Tc for the 3D EAL 
model in an equilibrium simulation, using a parallel tempering 
algorithm. '^''^^ To reach equilibrium between T = 1.0 and T = 
0.6, it was necessary to simulate 17 replicas of the system and 
a number of Monte Carlo sweeps of 2" with n = 2L + 8. At 
least up to 5 X 10^ samples for each lattice size were necessary 
to obtain a reliable disorder average. Due that the number 
of sweeps grow very fast with L, only lattices up to L = 8 
could be studied (probably this is due to the fact that Tc is 
lower than in the previous models). Figure |5] shows that the 
coiTelation length ^ /L (Ref. [191) and the Binder cumulant B 
(Ref. 123) cross at, respectively, ?;. = 0.8 1 ( 1 ) and ?;. = 0.79 (2) . 
Both values are compatible with the Tc obtained with our non- 
equilibrium method. 

The quantities estimates in this work for the three spin-glass 
models are shown in Table H] We report the values of expo- 
nent c (our best estimate of the true values of b and a), for 
which the error bar was determined taking the values of b and 
a as, respectively, the upper and the lower bounds of c. As we 
discuss before, the critical temperatures that we have calcu- 
lated here, agree very well with those obtained in equilibrium 
simulations. However, these Tc and the corresponding non- 
equilibrium critical exponents, differ from values reported in 
the literature: Tc = 1.19(1) and b = 0.056(3) for the EAB, 
Tc = 0.92(1) and b = 0.043(1) for the EAG, and Tc = 0.72(2) 
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FIG. 4: (Color online) Non-equilibrium exponents as function of 
temperature for (a) the 3D EAG and (b) the 3D EAL models, ob- 
tained for different ranges of t„ as indicated. 
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FIG. 5: (Color online) (a) Correlation length and (b) Binder cumulant 
as function of T for the 3D EAL model. Insets: data collapsing. 



TABLE I: Quantities calculated in this work. 



Parameter EAB 



EAG 



EAL 



Tc 1.135(5) 0.95(1) 0.815(5) 

c 0.044(2) 0.0455(5) 0.047(5) 

Xc/zc 0.29(1) 0.29(1) 0.254(2) 

Xo, 0.09(4) 0.09(1) 0.04(1) 

and b — 0.032(2) for the EAL model^ It is important to point 
out that at these same temperatures we obtain similar values 
ofb:b = 0.057, b = 0.0405, and b = 0.034 for, respectively, 
the EAB, the EAG and the EAL models. However, contrary to 
previous reports, we do not identify these temperatures with 
the Tc of each system. 

In addition, the values of exponent c that we have obtained 
are very similar to each other, and they are also quite close 
to b Ki 0.0464, the calculated value from the relation b = 
{d — 2 + r\)/2zc?^ where (i = 3 and we have used r\ ~ —0.375 
(Ref. 16) and Zc — 6.74.^^ As it has been shown previously 
in equilibrium simulations,^ these new results suggests that 
the universality class of the 3D Edwards-Anderson spin-glass 



model does not depend on the exact form of the bond dis- 
tribution. We can even notice that the coiTelation-Iength ex- 
ponents for the 3D EAL model, V = 2.7(2) and v = 2.5(2) 
that we have calculated from, respectively, the data collapse 
of the coiTelation length and the Binder cumulant (see insets 
in Fig.|5]), agree very well with those obtained previously for 
the 3D EAB and the 3D EAG models.^'^*^ 

On the other hand, the assumption of that the universality 
is violated, it is also based on quantities such as )ic/zc and 
Xoo = Hm,,,,^ooHm(^oop(r,fH.)/C(f ,f„.), the critical fluctuation- 
dissipation ratio. ^^ Our simulations show that for the 3D EAB 
and the 3D EAG models, the calculated values of Xc/zc andX„o 
agree, but differs from the coiTesponding one for the 3D EAL 
model (see Table |l]l. Although both quantities are believed 
to take universal values^ recently this conjecture has been 
questioned, showing that models belonging to the same uni- 
versality class at equilibrium, have different values of A^/zc 
orXoo.^^ To determine if these results are evidence of the non- 
universal character of Xc/zc and X„, or of the existence of 
different non-equilibrium universality classes, further investi- 
gations are required. ^^ 

In summary, we have proposed a non-equilibrium method 
to determine the Tc and the scaling exponents of the corre- 
lation and the integrated response functions. We apply this 
technique to the 3D Edwards-Anderson spin-glass model with 
three different bond distributions. The values of Tc that we 
obtain agree with those calculated in the equilibrium simula- 
tions. As the values of the exponent c are very close to each 
other, we conclude that for the 3D spin glasses the most prob- 
able scenario is that the universality class does not depend on 
the exact form of bond distribution. 
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